Projekat - Medicinski podaci NHANES

U ovom projektu korisceni su NHANES medicinski podaci. Posmatrani su podaci koji opisuju zivotne navike ljudi i njihovo zdravstveno stanje.

U prvom delu su pripremani podaci za statisticku analizu. Kasnije su ispitivane zavisnosti i potencijalne raspodele nekih obelezja. Na kraju su razmatrane neke od tehnika rukovanja nedostajucim vredostima u podacima.

library(foreign)
library(dplyr)
library(ggplot2)
library(mice)
library(EnvStats)
library(png)
library(grid)

#Dietary Interview - Individual Foods, First Day (DR1IFF_J)
DR1IFF_J <- read.xport("~/Desktop/Statistika/projekat/podaci/DR1IFF_J.XPT")
#Dietary Interview - Individual Foods, Second Day (DR2IFF_J)
DR2IFF_J <- read.xport("~/Desktop/Statistika/projekat/podaci/DR2IFF_J.XPT")

svi.prvi <- unique(DR1IFF_J[,"SEQN"])       #svi koji su radili anketu prvog dana
svi.drugi <- unique(DR2IFF_J[,"SEQN"])      #svi koji su radili anketu drugog dana

dorucak1 = filter(DR1IFF_J, DR1_030Z == 1 | DR1_030Z == 10) #1 i 10 su kodovi za doruckove na engleskom i spanskom
dorucak2 =filter(DR2IFF_J, DR2_030Z == 1 | DR2_030Z == 10)

unikatni1 = unique(dorucak1[,"SEQN"])     #svi koji su doruckovali prvog dana
unikatni2 = unique(dorucak2[,"SEQN"])     #svi koji su doruckovali drugog dana

oba.dana <- intersect(unikatni1,unikatni2)  #svi koji su radili anketu oba dana i doruckovali oba dana
prvi.dan <- setdiff(unikatni1, svi.drugi)   #svi koji su radili anketu samo prvi dan i doruckovali prvi dan

redovno.dorukcuju <- union(oba.dana,prvi.dan)
ne.doruckuju.redovno <- setdiff(union(svi.prvi,svi.drugi),redovno.dorukcuju) #nema missing values, pa moze razlika

redovno.dorukcuju <- data.frame(cbind(redovno.dorukcuju,T))
ne.doruckuju.redovno <- data.frame(cbind(ne.doruckuju.redovno,F))
names(redovno.dorukcuju) <- c("SEQN","redovno_doruckuje")
names(ne.doruckuju.redovno) <- c("SEQN","redovno_doruckuje")

podaci.dorucak <- rbind(redovno.dorukcuju,ne.doruckuju.redovno)
podaci.dorucak <- arrange(podaci.dorucak, SEQN)

######################
#Physical Activity (PAQ_J)

# PAQ605 - Vigorous work activity
# PAQ620 - Moderate work activity
# PAQ635 - Walk or bicycle
# PAQ650 - Vigorous recreational activities
# PAQ665 - Moderate recreational activities

PAQ_J <- read.xport("~/Desktop/Statistika/projekat/podaci/PAQ_J.XPT")
aktivni <- filter(PAQ_J, PAQ605 == 1 | PAQ620 == 1 | PAQ635 == 1 | PAQ650 == 1 | PAQ665 == 1) # 1 znaci odgovor je 'YES'
neaktivni <- setdiff(PAQ_J,aktivni)

aktivni <- data.frame(cbind(aktivni[,"SEQN"],TRUE))
neaktivni <- data.frame(cbind(neaktivni[,"SEQN"],FALSE))

names(aktivni) <- c("SEQN","aktivan")
names(neaktivni) <- c("SEQN","aktivan")

podaci.aktivni <- rbind(aktivni,neaktivni)
podaci.aktivni <- arrange(podaci.aktivni,SEQN)

#########################
#Body Measures (BMX_J)

# BMXBMI - Body Mass Index (kg/m**2)

BMX_J <- read.xport("~/Desktop/Statistika/projekat/podaci/BMX_J.XPT")
# funkcija proverava da li je bmi u normalnim granicama
proveri.bmi <- function(x){
  if(is.na(x)) return (x)
  if( x<=24.9 & x>=18.5)
    return (T) 
  else return (F)
}

podaci.bmi <- BMX_J[,c("SEQN","BMXBMI")]

podaci.bmi <- cbind(podaci.bmi,sapply(podaci.bmi[,2], proveri.bmi))
names(podaci.bmi)[2:3] <- c("bmi","normalan_bmi")

#########################
#Sleep Disorders (SLQ_J)

#SLD012 - Sleep hours - weekdays or workdays

SLQ_J <-read.xport("~/Desktop/Statistika/projekat/podaci/SLQ_J.XPT")

proveri.san <- function(x){
  if(is.na(x)) return (x)
  if( x<=8 & x>=7)
    return (T) 
  else return (F)
}

podaci.san <- SLQ_J[,c("SEQN","SLD012")]
podaci.san <- cbind(podaci.san,sapply(podaci.san[,2], proveri.san))
names(podaci.san)[2:3] <- c("san","optimalan_san")

#########################
#Alcohol Use (ALQ_J)
# ALQ111 - Ever had a drink of any kind of alcohol (2 znaci Yes)
# ALQ121 - Past 12 mo how often have alcohol drink
# ALQ130 - Avg # alcohol drinks/day - past 12 mos (prosecan broj pica za dane u kojima su pili alkohol)

ALQ_J <- read.xport("~/Desktop/Statistika/projekat/podaci/ALQ_J.XPT")
podaci.alkohol <- filter(ALQ_J, ALQ111==2 | ALQ121==0 | between(ALQ121,7,10))
ostali <- setdiff(ALQ_J,podaci.alkohol)

podaci.alkohol <- data.frame(cbind(podaci.alkohol[,"SEQN"],1))
names(podaci.alkohol) <- c("SEQN", "alkohol")

ostali <- data.frame(cbind(ostali[,"SEQN"],0))
names(ostali) <- c("SEQN", "alkohol")

podaci.alkohol <- rbind(podaci.alkohol,ostali)
podaci.alkohol<- arrange(podaci.alkohol, SEQN)


alkohol <- read.xport("~/Desktop/Statistika/projekat/podaci/ALQ_J.XPT")
alkohol <- alkohol[, c("SEQN","ALQ130","ALQ121","ALQ111")]
names(alkohol)[2:4] <- c("prosek_pica","koeficijent_redovnosti","ikad_pio")

obradi <- function(x){
  if(is.na(x[4])) return (NA)
  if(x[4] == 2) return (0)
  if(x[4] == 1){
    if(is.na(x[2])) return (NA)
    if(x[2]==777) return(NA)
    if(x[2]==999) return(NA)
    if(x[3]==0) return (0)
    if(x[3]==1) return (x[2]*365)
    if(x[3]==2) return (x[2]*365/1.5)
    if(x[3]==3) return (x[2]*365/2)
    if(x[3]==4) return (x[2]*365/3.5)
    if(x[3]==5) return (x[2]*52)
    if(x[3]==6) return (x[2]*12*2.5)
    if(x[3]==7) return (x[2]*12)
    if(x[3]==8) return (x[2]*9)
    if(x[3]==9) return (x[2]*4.5)
    if(x[3]==10) return (x[2]*1.5)
    return (NA)
  }
}

ukupno_pica <- apply(alkohol,1,obradi)
alkohol <- cbind(alkohol, ukupno_pica)
names(alkohol)[5] <- "ukupno_pica_godisnje"
#alkohol <- arrange(alkohol,desc(ukupno_pica_godisnje))
alkohol <- alkohol[,c("SEQN","ukupno_pica_godisnje")]

##########################
#Smoking - Cigarette Use (SMQ_J)

# SMQ020 - Smoked at least 100 cigarettes in life
# SMQ040 - Do you now smoke cigarettes?

SMQ_J <- read.xport("~/Desktop/Statistika/projekat/podaci/SMQ_J.XPT")
nepusaci <- filter(SMQ_J, SMQ020 ==2 | SMQ040 ==3)
pusaci <- filter(SMQ_J, SMQ040 ==1 | SMQ040 ==2)

nepusaci <- data.frame(cbind(nepusaci[,"SEQN"],T))
pusaci <- data.frame(cbind(pusaci[,"SEQN"],F))

names(nepusaci) <- c("SEQN","nepusac")
names(pusaci) <- c("SEQN","nepusac")

podaci.cigare <- rbind(nepusaci, pusaci)

############################
#Demographic Variables and Sample Weights (DEMO_J)
# RIAGENDR - Gender
# RIDAGEYR - Age in years at screening

demografija <- read.xport("~/Desktop/Statistika/projekat/podaci/DEMO_J.XPT")
demografija <- demografija[, c("SEQN","RIAGENDR", "RIDAGEYR")]
names(demografija)[2:3] <- c("pol","godine")

##############################
#Blood Pressure (BPX_J)

# BPXCHR - 60 sec HR (30 sec HR * 2)  (puls - deca do 7 godina)
# BPXPLS - 60 sec. pulse (30 sec. pulse * 2) (puls - stariji od 7 godina)
# BPXSY1 - Systolic: Blood pres (1st rdg) mm Hg
# BPXSY2 - Systolic: Blood pres (2nd rdg) mm Hg
# BPXDI1 - Diastolic: Blood pres (1st rdg) mm Hg
# BPXDI2 - Diastolic: Blood pres (2nd rdg) mm Hg

krvni_pritisak <- read.xport("~/Desktop/Statistika/projekat/podaci/BPX_J.XPT")
krvni_pritisak <- krvni_pritisak[,c("SEQN","BPXCHR", "BPXPLS","BPXSY1","BPXSY2","BPXDI1","BPXDI2")]
krvni_pritisak[,3] <- ifelse(is.na(krvni_pritisak$BPXPLS),krvni_pritisak$BPXCHR,krvni_pritisak$BPXPLS)
krvni_pritisak[,4] <- ifelse(is.na(krvni_pritisak$BPXSY1),krvni_pritisak$BPXSY2,krvni_pritisak$BPXSY1)
krvni_pritisak[,6] <- ifelse(is.na(krvni_pritisak$BPXDI1),krvni_pritisak$BPXDI2,krvni_pritisak$BPXDI1)
krvni_pritisak <- krvni_pritisak[,c(1,3,4,6)]
names(krvni_pritisak)[2:4] <- c("puls","gornji_pritisak","donji_pritisak")


###############################
#Cholesterol - Total (TCHOL_J)

#LBXTC - Total Cholesterol (mg/dL)

holesterol <- read.xport("~/Desktop/Statistika/projekat/podaci/TCHOL_J.XPT")
holesterol <- holesterol[,c("SEQN","LBXTC")]
names(holesterol)[2] <- "holesterol"

###############################
#Current Health Status (HSQ_J)

#HSD010 - General health condition

HSQ_J <- read.xport("~/Desktop/Statistika/projekat/podaci/HSQ_J.XPT")
podaci.zdravlje <- HSQ_J[,c("SEQN","HSD010")]
names(podaci.zdravlje)[2] <- "osecaj_zdravlja"
#namestamo da za skalu osecaja trenutnog zdravlja od 1 do 5, najbolji osecaj bude 5
podaci.zdravlje[,2] = sapply(podaci.zdravlje[,2],function(x) ifelse(is.na(x) | x==9 | x==7, NA, 6-x))

Spajamo sve podatke u jednu tabelu.

tabela <- merge(demografija, krvni_pritisak,by = "SEQN", all = T)
tabela <- merge(tabela, holesterol,by = "SEQN", all = T)
tabela <- merge(tabela, alkohol,by = "SEQN", all = T)
tabela <- merge(tabela, podaci.cigare,by = "SEQN", all = T)
tabela <- merge(tabela,podaci.alkohol, by = "SEQN", all = T)
tabela <- merge(tabela,podaci.san, by = "SEQN", all = T)
tabela <- merge(tabela,podaci.bmi, by = "SEQN", all = T)
tabela <- merge(tabela,podaci.aktivni, by = "SEQN", all = T)
tabela <- merge(tabela,podaci.dorucak, by = "SEQN", all = T)
tabela <- merge(tabela,podaci.zdravlje, by = "SEQN", all = T)

# prebacujemo vrednosti indikatorskih kolona u True i False zbog plotovanja
prebaci.u.T.F <- function(x){
  if(is.na(x)) return (NA)
  if(x==1) return (TRUE)
  return (FALSE)
}

prebaci.kolonu <- function(x){
  sapply(x, prebaci.u.T.F)
}

tabela[,c("nepusac","alkohol","optimalan_san","aktivan", "redovno_doruckuje","normalan_bmi")] = 
    apply(tabela[,c("nepusac","alkohol","optimalan_san","aktivan", "redovno_doruckuje","normalan_bmi")], 2, prebaci.kolonu)

head(tabela)

##    SEQN pol godine puls gornji_pritisak donji_pritisak holesterol
## 1 93703   2      2  120              NA             NA         NA
## 2 93704   1      2  114              NA             NA         NA
## 3 93705   2     66   52              NA             NA        157
## 4 93706   1     18   82             112             74        148
## 5 93707   1     13  100             128             38        189
## 6 93708   2     66   68             138             78        209
##   ukupno_pica_godisnje nepusac alkohol  san optimalan_san  bmi normalan_bmi
## 1                   NA      NA      NA   NA            NA 17.5        FALSE
## 2                   NA      NA      NA   NA            NA 15.7        FALSE
## 3                   12    TRUE    TRUE  8.0          TRUE 31.7        FALSE
## 4                    0    TRUE    TRUE 10.5         FALSE 21.5         TRUE
## 5                   NA      NA      NA   NA            NA 18.1        FALSE
## 6                    0    TRUE    TRUE  8.0          TRUE 23.7         TRUE
##   aktivan redovno_doruckuje osecaj_zdravlja
## 1      NA                NA              NA
## 2      NA              TRUE              NA
## 3    TRUE              TRUE               3
## 4    TRUE             FALSE               4
## 5      NA              TRUE               3
## 6    TRUE              TRUE               3

Posmatramo odgovor ispitanika kako se oseca fizicki, tj. da li se oseca zdavo, u odnosu na broj zdravih navika koje upraznjava. Zdrave navike koje posmatramo jesu: Da li osoba ne pusi, ne pije vise od jedanput mesecno, bar umereno je fizicki aktivna, ima optimalan san (7-8 sati), da li redovno doruckuje i da li ima bmi index u normalnim granicama.

broj_upraznjavanja_navika <- tabela$nepusac+tabela$alkohol+tabela$aktivan+tabela$optimalan_san+tabela$redovno_doruckuje+tabela$normalan_bmi
tabela <- cbind(tabela,broj_upraznjavanja_navika)


ggplot(na.omit(tabela[,c("broj_upraznjavanja_navika","osecaj_zdravlja")]), aes(x=factor(broj_upraznjavanja_navika), fill=factor(osecaj_zdravlja)))+
  geom_bar(stat="count", width=0.6) + 
  geom_text(stat = "count", aes(label=..count..), size = 3, position = position_stack(0.6))+
  theme_minimal()

by_br_navika <- group_by(na.omit(tabela[,c("broj_upraznjavanja_navika","osecaj_zdravlja")]),broj_upraznjavanja_navika)
by_br_navika <- summarise(by_br_navika, osecaj_zdravlja_mean = mean(osecaj_zdravlja))
rezultat_navika2 <- data.frame(by_br_navika)
by_br_navika <- group_by(na.omit(tabela[,c("broj_upraznjavanja_navika","osecaj_zdravlja")]),broj_upraznjavanja_navika)
by_br_navika <- summarise(by_br_navika,osecaj_zdravlja_sd = sd(osecaj_zdravlja))
rezultat_navika2 <- cbind(rezultat_navika2, data.frame(by_br_navika[2]))
rezultat_navika2

##   broj_upraznjavanja_navika osecaj_zdravlja_mean osecaj_zdravlja_sd
## 1                         0             2.625000          1.1877349
## 2                         1             2.870504          0.9075731
## 3                         2             2.896030          0.9626369
## 4                         3             3.070463          0.9428182
## 5                         4             3.183134          0.9534862
## 6                         5             3.358011          0.9205162
## 7                         6             3.622093          0.9983497

tabela[,"alkohol"] = sapply(tabela[,"alkohol"], function(x) ifelse(x,"ne pije","pije"))
tabela[,"pol"] = sapply(tabela[,"pol"], function(x) ifelse(x==1,"Musko","Zensko"))
tabela[,"nepusac"] = sapply(tabela[,"nepusac"], function(x) ifelse(x,"ne pusi","pusi"))

Procenat ljudi koji piju (cesce od jednom mesecno) u odnosu na godine (18god - 80+god). Vidimo da najvise ljudi piju u dobi od 22. godine do 30. godine (vise od 50% ljudi pije).

ggplot(na.omit(tabela[,c("godine","alkohol")]), aes(x=factor(godine), fill=factor(alkohol)))+
  geom_bar(stat="count", width=0.6, position = "fill") +
  labs(title="Procenat ljudi koji piju u odnosu na godine", x= "18 - 80+ godina", y = "procenat ljudi" )

ggplot(na.omit(tabela[,c("pol","nepusac")]), aes(x=pol, fill=nepusac)) + geom_bar()+ 
  geom_text(stat = "count", aes(label=..count..), size = 5, position = position_stack(0.5))+
  labs(title="Broj pusaca u odnosu na pol")

ggplot(na.omit(tabela[,c("nepusac","alkohol")]), aes(x=nepusac, fill=alkohol)) +
  geom_bar(stat="count", width=0.6, position = "fill")+
  labs(title="Zavisnost da li ljudi piju od toga da li puse")

ggplot(na.omit(tabela[,c("SEQN","alkohol","nepusac")]), aes(x=alkohol, fill=nepusac)) +
  geom_bar(stat="count", width=0.6, position = "fill") +
  labs(title="Zavisnost da li ljudi puse od toga da li piju")

ggplot(na.omit(tabela[,c("SEQN","nepusac","redovno_doruckuje")]), aes(x=nepusac, fill=redovno_doruckuje)) +
  geom_bar(stat="count", width=0.6, position = position_dodge()) +
  geom_text(stat = "count", aes(label=..count..), size = 5, position =position_dodge(0.6) ,vjust = 1.6)+
  labs(title="Zavisnost da li ljudi doruckuju od toga da li puse")

ggplot(na.omit(tabela[,c("pol","nepusac")]), aes(x=pol, fill=nepusac)) + geom_bar()+ 
  geom_text(stat = "count", aes(label=..count..), size = 5, position = position_stack(0.5))+
  labs(title="Broj pusaca u odnosu na pol")

Ilustracija odnosa broja ljudi koji piju u odnosu na pol. Radimo Hi-kvadrat test nezavnosti za ova 2 obelezja.

ggplot(na.omit(tabela[,c("pol","alkohol")]), aes(x=pol, fill=alkohol)) + geom_bar(stat="count", width=0.6, position = "fill") + 
  labs(title="Procenat ljudi koji piju vise od jedanput mesecno u odnosu na pol")

podaci.test.nez <- na.omit(tabela[,c("pol","alkohol")])
l <- nrow(podaci.test.nez)
podaci.test.nez <-c(nrow(filter(podaci.test.nez, pol == "Musko" & alkohol == "ne pije")),
                    nrow(filter(podaci.test.nez, pol == "Musko" & alkohol == "pije")),
                    nrow(filter(podaci.test.nez, pol == "Zensko" & alkohol == "ne pije")),
                    nrow(filter(podaci.test.nez, pol == "Zensko" & alkohol == "pije")))
podaci.test.nez <- matrix(podaci.test.nez, ncol = 2, byrow = T)
p_redovi <- rowSums(podaci.test.nez)/l
p_kolone <- colSums(podaci.test.nez)/l
vrednosti <- expand.grid(p_redovi,p_kolone)
vrednosti <- vrednosti[1] * vrednosti[2]
vrednosti <- matrix(vrednosti[[1]], ncol = 2)* l

test_statistika <- sum((podaci.test.nez - vrednosti)^2 / vrednosti)

test_statistika ima chisq raspodelu sa (2-1)*(2-1)= 1 stepenom slobode. Proveravamo da li je statistika upala u kriticnu oblast i posto jeste odbacujemo hipotezu da su obelezja nezavisna.

test_statistika > qchisq(0.99,1)

## [1] TRUE

Pokusacemo da odredimo raspodelu za puls kod odraslih ljudi

pulsevi <-filter(tabela,godine>=18 & !is.na(puls))$puls
#ocenjujemo nepoznate parametre metodom maksimalne verodostojnosti
egamma(pulsevi)

## 
## Results of Distribution Parameter Estimation
## --------------------------------------------
## 
## Assumed Distribution:            Gamma
## 
## Estimated Parameter(s):          shape = 37.501935
##                                  scale =  1.916804
## 
## Estimation Method:               MLE
## 
## Data:                            pulsevi
## 
## Sample Size:                     5269

# ocenjeni parametri
# shape     scale 
# 37.501935  1.916804 

#delimo uzorak na kategorije i racunamo broj elemanata koji upadaju u njih
n<- length(pulsevi)
k = floor(log2(n))+1
kvant <- quantile(pulsevi, probs = seq(0,1,1/k), na.rm = T)
mk = rep(0,k)

kvant[length(kvant)] = Inf
kvant[1] = 0

for( i in 1:k){
  mk[i] = length(pulsevi[kvant[i]<pulsevi & pulsevi<=kvant[i+1]])
}

# Trazimo parametre shape i scale td. funkcija XkvadratNorm ima minimalnu vrednost
# Za to koristimo funkciju 'optim' koja koristi numericki pristup za ovaj problem
XkvadratNorm <- function(param, a, M, pr=F)
{
  oblik=param[1]
  skala=param[2]
  k=length(a)
  n=sum(M)
  E=(pgamma(a[2:k] ,oblik, scale = skala)-pgamma(a[1:(k-1)],oblik, scale = skala))*n
  X2=sum((M-E)^2/E)
  if (pr) print (E)
  return (X2)
}
param = c(37.501935,1.916804)
param1 = optim(param, fn=XkvadratNorm, a=kvant, M=mk)$par
#dobijeni parametri su
param1

## [1] 37.442845  1.884889

hist(pulsevi,probability = T, ylim = c(0,0.04),main = "Histogram pulsa odraslog coveka")
curve(dgamma(x,37.424135,scale = 1.885739),add = T,col = "darkgreen", lwd = 2)
curve(dnorm(x,mean = mean(pulsevi), sd = sd(pulsevi)), col = "brown",add = T,lwd = 2)
legend("topright", legend = c("Gama", "Normalna"), fill = c("darkgreen","red"))

hikv = XkvadratNorm(param1,kvant,mk)
#hikv = 63.03893; qgamma(0.99,k-3) = 18.78312
# k broj kategorija i 2 ocenjena parametra
#statistika upala u kriticnu oblast => nije gama raspodela
hikv > qgamma(0.99,k-3)

## [1] TRUE

=>Odbacujemo hipotezu da puls kod odraslih ljudi ima Gama raspodelu

Sledece obelezje za koje proveravamo raspodelu je broj alkoholnih pica koje osoba popije u toku godine.

pica <- filter(tabela,!is.na(ukupno_pica_godisnje))$ukupno_pica_godisnje

epareto(pica+1, method = "mle")

## 
## Results of Distribution Parameter Estimation
## --------------------------------------------
## 
## Assumed Distribution:            Pareto
## 
## Estimated Parameter(s):          location = 1.0000000
##                                  shape    = 0.3091797
## 
## Estimation Method:               mle
## 
## Data:                            pica + 1
## 
## Sample Size:                     4071

# $parameters
# location     shape 
# 1.0000000 0.3091797 
hist(pica[pica<1500]+1, probability = T,main = "Histogram popijenih pica u roku od godinu dana\n(Prikazano do 1500)", xlab = "popijeno pica", ylab = "dpareto") #krece od 0
curve(dpareto(x,1,0.3091797),add = T)

alfa = 0.3091797

k=6
kvant <- quantile(pica, probs = seq(0,1,1/k), na.rm = T)
Npk = rep(0,k)
Mk = rep(0,k)
n<- length(pica)
for( i in 1:k){
  Npk[i] = (ifelse(i==k,1,ppareto(kvant[i+1]+1,1,alfa)) - ppareto(kvant[i]+1,1,alfa))*n
}
Mk[1] = length(pica[pica==min(pica)])
for( i in 1:k){
  Mk[i] = Mk[i]+length(pica[kvant[i]<pica & pica<=kvant[i+1]])
}
rbind(Mk,Npk)

##         [,1]     [,2]     [,3]     [,4]     [,5]   [,6]
## Mk   909.000  669.000 475.0000 761.0000 584.0000 673.00
## Npk 1004.336 1069.004 492.8327 495.5614 280.6052 728.66

Hikv = 0

for (i in 1:5){
  Hikv = Hikv + (Mk[i]-Npk[i])^2/Npk[i]
}
#Hikv = 629.5832; qchisq=13.2767
Hikv > qchisq(0.99,k-2)

## [1] TRUE

Za gornji pritisak radimo kao i za puls i dobijamo isti ishod testiranja hipoteze.

gor_prit <- filter(tabela,!is.na(gornji_pritisak))$gornji_pritisak
egamma(gor_prit)

## 
## Results of Distribution Parameter Estimation
## --------------------------------------------
## 
## Assumed Distribution:            Gamma
## 
## Estimated Parameter(s):          shape = 38.066622
##                                  scale =  3.202626
## 
## Estimation Method:               MLE
## 
## Data:                            gor_prit
## 
## Sample Size:                     6678

n<- length(gor_prit)
k = floor(log2(n))+1
kvant <- quantile(gor_prit, probs = seq(0,1,1/k), na.rm = T)
mk = rep(0,k)

kvant[length(kvant)] = Inf
kvant[1] = 0
mk[1] = length(gor_prit[gor_prit==min(gor_prit)])
for( i in 1:k){
  mk[i] = mk[i]+length(gor_prit[kvant[i]<gor_prit & gor_prit<=kvant[i+1]])
}

XkvadratNorm <- function(param, a, M, pr=F)
{
  oblik=param[1]
  skala=param[2]
  k=length(a)
  n=sum(M)
  E=(pgamma(a[2:k] ,oblik, scale = skala)-pgamma(a[1:(k-1)],oblik, scale = skala))*n
  X2=sum((M-E)^2/E)
  if (pr) print (E)
  return (X2)
}
param = c(38.066622,3.202626)
param1 = optim(param, fn=XkvadratNorm, a=kvant, M=mk)$par

## Warning in pgamma(a[2:k], oblik, scale = skala): NaNs produced

## Warning in pgamma(a[1:(k - 1)], oblik, scale = skala): NaNs produced

param1

## [1] 37.137148  3.234645

hist(gor_prit,probability = T, ylim = c(0,0.04), main = "Histogram gornjeg pritiska")
curve(dgamma(x,param1[1],scale = param[2]),add = T,col = "darkgreen", lwd = 2)
curve(dnorm(x,mean = mean(gor_prit), sd = sd(gor_prit)), col = "brown",add = T,lwd = 2)
legend("topright", legend = c("Gama", "Normalna"), fill = c("darkgreen","red"))

hikv = XkvadratNorm(param1,kvant,mk, pr = F)
#hikv = 345.7659; qgamma(0.99,k-3) = 18.78312
hikv > qgamma(0.99,k-3)

## [1] TRUE

Nedostajuce vrednosti u bazi

Podaci koji nedostaju mogu predstavljati ne tako trivijalan problem kada se analizira skup podataka I upravljanje njima obično nije tako jednostavno. Ako je količina podataka koji nedostaju vrlo mala u odnosu na veličinu skupa podataka, izostavljanje nekoliko uzoraka sa nedostajućim karakteristikama može biti najbolja strategija kako ne bi došlo do pristrasnosti u analizi, međutim izostavljanje dostupnih tačaka podataka uskraćuje informacije o podacima. U zavisnosti od situacije sa kojom smo suočeni, možda je bolje da potražimo druge nacine pre nego što izbrišemo potencijalno korisne podatke iz skupa podataka.

—Postoje tri vrste podataka koji nedostaju:

Posmatrajmo nedostajuce podatke u nekom atributu Y. Mozemo oznaciti sa Ya podskup Y onih podataka koji su dostupni, I sa Ym podskup Y onih podataka koji su nedostajuci. Ostale podatke iz drugih atributa oznacicemo sa X. M je indikator da li neki podatak iz Y nedostaje.

MCAR: missing completely at random. Ovo je poželjan scenario u slučaju da nedostaju podaci. Činjenica da nedostaje određena vrednost nema nikakve veze sa njenom hipotetičkom vrednošću i vrednostima drugih promenljivih. P(M| Y,X) = P(M)

MAR: missing at random. Podaci slucajno nedostaju, odnosno to sto podaci nedostaju nije povezano sa samim podacima koji nedostaju, ali je povezano sa nekim od posmatranih podataka. Odnosno vazi: P(M| Y,X) = P(M| Ya,X)

MNAR: missing not at random. Nedostatak neslučajnih podataka je ozbiljniji problem i u ovom slučaju bi bilo pametno dodatno proveriti postupak prikupljanja podataka i pokušati shvatiti zašto informacije nedostaju. Na primer, ako većina ljudi u anketi nije odgovorila na određeno pitanje, zašto su to učinili? Da li je pitanje bilo nejasno? Dva moguća razloga su da vrednost koja nedostaje zavisi od hipoteticke vrednosti (npr. ljudi sa visokim platama uglavnom ne žele da otkriju svoje prihode u anketama) ili da vrednost koja nedostaje zavisi od vrednosti neke druge promenljive (npr. pretpostavimo da žene uglavnom ne žele da otkriju njihove godine). Ovde na nedostajuću promenljivu starosti utiče vrednost promenljive pola. P(M| Y,X) = P(M| Ym ,Ya,X)

Tehnike postupanja sa nedostajucim podacima:

Podaci koji nedostaju smanjuju reprezentativnost uzorka i stoga mogu iskriviti zaključke o populaciji. Uopšteno govoreći, postoje tri glavna pristupa za rukovanje podacima koji nedostaju: (1) brisanje - gde se uzorci sa nedostajucim podacima odbacuju iz dalje analize, (2) imputacija - gde se vrednosti popunjavaju na mestu nedostajućih podataka i (3) analiza - direktnom primenom metoda na koje vrednosti koje nedostaju ne utiču.

Metode koje redukuju skupove podataka tako sto izbacuju nedostajuce vrednosti su:

Listwise deletion (Brisanje svih elemenata iz uzorka kojima fali bar jedna vrednost) • Podrazumevani metod za vecinu statističkih paketa. • Vrši se za slučaj da podaci neodostaju po tipu MCAR i kada je procenat nedostajucih podataka mali. • Nedostatak: Dimenzije posmatranog skupa se mogu značajno smanjiti • Primenjivati kada je nedostajucih podataka <5%
Pairwise deletion (Brisanje u paru, npr. kada poredimo dva atributa koristimo sve elemente koje poseduju oba podatka bez obzira na postojanje podataka iz ostalih atributa) • Isključivanje ili brisanje samo nedostajucih podataka • Isključiti iz analize samo opservaciju za koju nedostaje vrednost konkretne varijable
Deleting Columns (Brisanje kolona) • Isključiti iz analize varijablu za koju nedostaje visok procenat podataka
• Nedostatak: može se isključiti varijabla koja je od značaja za dalju analizu

Imputation - Popunjavanje nedostajucih podataka

U statistici je imputacija postupak zamene nedostajućih podataka nekim vrednostima. Postoje tri glavna problema koja nedostajuci podaci uzrokuju: podaci koji nedostaju mogu uvesti značajnu količinu pristrasnosti, učiniti rukovanje i analizu podataka složenijim i stvoriti smanjenje efikasnosti. Budući da podaci koji nedostaju mogu stvoriti probleme za analizu podataka, imputiranje se smatra načinom da se izbegnu nedostaci koje proizvodi brisanje slučajeva koji imaju vrednosti koje nedostaju. Odnosno, kada nedostaje jedna ili više vrednosti za slucaj, većina statističkih paketa odbacuje svaki slucaj koji ima vrednost koja nedostaje, što može uvesti pristrasnost ili uticati na reprezentativnost rezultata. Imputacija čuva sve slučajeve zamenom nedostajućih podataka procenjenom vrednošću na osnovu drugih dostupnih informacija. Jednom kada su imputirane sve vrednosti koje nedostaju, skup podataka može se analizirati pomoću standardnih tehnika za kompletne podatke. Neke od poznatih tehnika za rešavanje problema podataka koji nedostaju su: hot deck i cold deck imputacija, imputacija srednje vrednosti; nenegativna matricna faktorizacija; regresiona imputacija, stohasticka imputacija i visetruka imputacija.

–Višestruke imputacije:

Da bi se prevazisao problem velike disperzije, Rubin(1987) je razvio metod za uzimanje prosecnih vrednosti iz vise ponovljenih imputacija skupa podataka. Nedostajuci podaci se popunjavaju kroz slučajni proces koji uključuje neizvesnost. Popunjavanje se vrši N puta, u cilju kreiranje N “potpunih” skupova podataka. Parametri se procenjuju svaki put uključujuci i standardnu grešku..Prosečna vrednost svih N procena se kombinuje da bi se dobila nedostajuca vrednost.

Sve metode višestruke imputacije se sastoje iz sledeca tri koraka:

Imputacija - Slično kao kod pojedinačne imputacije, vrednosti koje nedostaju se pripisuju. Međutim, imputirane vrednosti se crpe m puta iz distribucije, a ne samo jednom. Na kraju ovog koraka trebalo bi da postoji m kompletiranih skupova podataka.
Analiza - Analizira se svaki od m skupova podataka. Na kraju ovog koraka trebalo bi da postoji m analiza.
Udruživanje - m rezultati su objedinjeni u jedan rezultat izračunavanjem srednje vrednosti, varijanse i intervala poverenja za svaku vrednost ili kombinovanjem simulacija iz svakog odvojenog modela.

Kao što postoji više metoda pojedinačne imputacije, tako postoji i više metoda višestruke imputacije. Jedna prednost koju višestruka imputacija ima u odnosu na pojedinačnu imputaciju i kompletne metode slučaja je ta što je višestruka imputacija fleksibilna i može se koristiti u širokom spektru scenarija. Višestruka imputacija može se koristiti u slučajevima kada podaci nedostaju potpuno nasumično (MCAR), nedostaju nasumično (MAR), pa čak i kada podaci neslucajno nedostaju (MNAR). Međutim, primarni metod višestruke imputacije je višestruka imputacija lančanim jednačinama (multiple imputation by chained equations - MICE). Takođe je poznato kao „potpuno uslovna specifikacija“ i, „višestruka imputacija sekvencijalne regresije“. Pokazalo se da MICE vrlo dobro radi kod nedostajanja nasumičnih podataka, mada postoje dokazi koji kroz simulacionu studiju sugerišu da sa dovoljnim brojem pomoćnih promenljivih može raditi i na podacima koji neslucajno nedostaju.

Kao što je navedeno u prethodnom odeljku, pojedinačna imputacija ne uzima u obzir gresku u imputovanim vrednostima. Nakon imputacije, podaci se tretiraju kao da su stvarne vrednosti. Zanemarivanje greske u imputaciji može i dovešće do grešaka u zaključcima. Imputiranjem više puta, višestruka imputacija smanjuje nesigurnost i opseg vrednosti koje je prava vrednost mogla uzeti.

Iako je pojedinačnu imputaciju lakše primeniti, ni višestruku imputaciju nije previse teško primeniti, a pokazuje se kao efikasniji nacin. Postoji širok spektar različitih statističkih paketa u različitim statističkim softverima koji lako omogućavaju nekome da izvrši višestruku imputaciju. Na primer, MICE paket omogućava korisnicima u R-u da izvrše višestruku imputaciju pomoću MICE metode.

*Paket mice pruža funkciju md.pattern() za bolje razumevanje obrasca podataka koji nedostaju:

Ona nam pokazuje vise informacija o tome koji podataci fale, sto nam moze koristiti za dalju analizu podataka i odabir metoda za rukovanje njima.

library(mice)
md.pattern(tabela[,c("alkohol","aktivan","nepusac","normalan_bmi","optimalan_san","redovno_doruckuje")],rotate.names = T)

##      normalan_bmi redovno_doruckuje optimalan_san aktivan nepusac alkohol      
## 4889            1                 1             1       1       1       1     0
## 280             1                 1             1       0       0       0     3
## 35              1                 1             0       1       1       1     1
## 1837            1                 1             0       0       0       0     4
## 506             1                 0             1       1       1       1     1
## 12              1                 0             1       0       0       0     4
## 4               1                 0             0       1       1       1     2
## 442             1                 0             0       0       0       0     5
## 56              0                 1             1       1       1       1     1
## 3               0                 1             0       1       1       1     2
## 541             0                 1             0       0       0       0     5
## 38              0                 0             1       1       1       1     2
## 320             0                 0             1       1       1       0     3
## 12              0                 0             1       0       0       0     5
## 2               0                 0             0       1       1       1     3
## 3               0                 0             0       1       1       0     4
## 274             0                 0             0       0       0       0     6
##              1249              1613          3141    3398    3398    3721 16520

Primenicemo metod MICE na neke atribute iz pocetne tabele, kako bi resili problem NA vrednosti:

nova_tabela = tabela[,c("godine","puls","gornji_pritisak","donji_pritisak","holesterol")]
#md.pattern(nova_tabela)

head(nova_tabela)

##   godine puls gornji_pritisak donji_pritisak holesterol
## 1      2  120              NA             NA         NA
## 2      2  114              NA             NA         NA
## 3     66   52              NA             NA        157
## 4     18   82             112             74        148
## 5     13  100             128             38        189
## 6     66   68             138             78        209

sapply(nova_tabela, function(x) sum(is.na(x))) # proveravamo koliko ima NA vrednosti za svaki atribut

##          godine            puls gornji_pritisak  donji_pritisak      holesterol 
##               0             973            2576            2576            2516

str(nova_tabela) # prikaz strukture skupa podataka

## 'data.frame':    9254 obs. of  5 variables:
##  $ godine         : num  2 2 66 18 13 66 75 0 56 18 ...
##  $ puls           : num  120 114 52 82 100 68 74 146 62 68 ...
##  $ gornji_pritisak: num  NA NA NA 112 128 138 120 NA 108 112 ...
##  $ donji_pritisak : num  NA NA NA 74 38 78 66 NA 68 68 ...
##  $ holesterol     : num  NA NA 157 148 189 209 176 NA 238 182 ...

inputed <- mice(nova_tabela,m=5,maxit=50,meth='pmm',seed = 100,print = FALSE) # predictive mean matching metod za imputaciju

inputed$predictorMatrix # matrica predikcije

##                 godine puls gornji_pritisak donji_pritisak holesterol
## godine               0    1               1              1          1
## puls                 1    0               1              1          1
## gornji_pritisak      1    1               0              1          1
## donji_pritisak       1    1               1              0          1
## holesterol           1    1               1              1          0

inputed$method # metode za imputaciju nedostajucih vrednosti za svaki atribut

##          godine            puls gornji_pritisak  donji_pritisak      holesterol 
##              ""           "pmm"           "pmm"           "pmm"           "pmm"

head(inputed$imp$holesterol) # provera imputiranih podataka za odredjen atribut npr. holesterol

##      1   2   3   4   5
## 1  151 172 216 149 164
## 2  138 199 128 148 170
## 8  161 136 199 140 197
## 18 159 131 122 153 133
## 22 173 132 129 113 165
## 23 153 211 152 146 144

summary(inputed)

## Class: mids
## Number of multiple imputations:  5 
## Imputation methods:
##          godine            puls gornji_pritisak  donji_pritisak      holesterol 
##              ""           "pmm"           "pmm"           "pmm"           "pmm" 
## PredictorMatrix:
##                 godine puls gornji_pritisak donji_pritisak holesterol
## godine               0    1               1              1          1
## puls                 1    0               1              1          1
## gornji_pritisak      1    1               0              1          1
## donji_pritisak       1    1               1              0          1
## holesterol           1    1               1              1          0

completedData <- complete(inputed,4)
sapply(completedData, function(x) sum(is.na(x)))

##          godine            puls gornji_pritisak  donji_pritisak      holesterol 
##               0               0               0               0               0

densityplot(inputed, main="Uporedni prikaz gustina posmatranih i imputiranih podataka")

# pregled svih metoda za imputaciju u R-u:
methods(mice)

## Warning in .S3methods(generic.function, class, envir): function 'mice' appears
## not to be S3 generic; found functions that look like S3 methods

##  [1] mice.impute.2l.bin       mice.impute.2l.lmer      mice.impute.2l.norm     
##  [4] mice.impute.2l.pan       mice.impute.2lonly.mean  mice.impute.2lonly.norm 
##  [7] mice.impute.2lonly.pmm   mice.impute.cart         mice.impute.jomoImpute  
## [10] mice.impute.lda          mice.impute.logreg       mice.impute.logreg.boot 
## [13] mice.impute.mean         mice.impute.midastouch   mice.impute.mnar.logreg 
## [16] mice.impute.mnar.norm    mice.impute.norm         mice.impute.norm.boot   
## [19] mice.impute.norm.nob     mice.impute.norm.predict mice.impute.panImpute   
## [22] mice.impute.passive      mice.impute.pmm          mice.impute.polr        
## [25] mice.impute.polyreg      mice.impute.quadratic    mice.impute.rf          
## [28] mice.impute.ri           mice.impute.sample       mice.mids               
## [31] mice.theme              
## see '?methods' for accessing help and source code

Na sledecoj slici imamo ilustraciju kako radi mice algoritam u toku jedne iteracije. Na pocetku inicijalizujemo nedostajuce vrednostima, sa npr. srednjom vrednoscu ili nekom slucajnom velicinom. Posle menjamo inicijalno nedostajucu vrednost posevno nekom metodom (npr. linearnom regeresijom), pritom koristeci ostale vrednosti i koje su dodate. Isto uraditi i za sve ostale inicijalno nedostajuce vrednosti iduci nekim redom, npr. s leva na desno.

Iterativno se izvrsavaju ovi koraci, dok se vrednosti koje zelimo da dodamo dovoljno ne stabilizuju ili dok ne izvrsimo neki broj iteracija.

Projekat - Medicinski podaci NHANES

Statistika za automatsku analizu podataka - Industrija 4.0, Matematički fakultet i Mašinski fakultet, Univerzitet u Beogradu

Stefan Kovač & Katarina Pantović

12.1.2021

Nedostajuce vrednosti u bazi

Tehnike postupanja sa nedostajucim podacima:

Imputation - Popunjavanje nedostajucih podataka